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Abstract:  A  depth-averaged  2-D  model  has  been  developed  for  simulating 
coastal  barrier  breaching  processes.  The  model  solves  the  generalized  shallow 
water  equations  and  the  non-equilibrium  sediment  transport  and  bed  change 
equations  in  a  coupled  fashion  using  an  explicit  finite  volume  method  on  a 
rectangular  grid.  It  considers  interactions  among  transient  flow,  strong  sediment 
transport,  and  rapid  bed  change  by  including  bed  change  and  variable  flow  density 
in  the  flow  continuity  and  momentum  equations.  The  model  adopts  the  HLL 
approximate  Riemann  solver  to  handle  the  mixed-regime  flows  near  the  breach 
and  the  wetting-drying  problem.  The  developed  model  has  been  tested  against  a 
lab  experiment  of  dike-break  flow  over  movable  bed  and  a  field  experiment  of  sea 
dike  breaching.  The  simulation  results  demonstrate  that  the  model  is  capable  of 
calculating  the  initial  development  stage  of  coastal  barrier  breaching  under 
supercritical  flow.  Future  work  will  improve  the  model  to  include  waves  for  more 
general  application  in  the  coastal  context. 

Introduction 

Breaching  often  occurs  on  coastal  barrier  islands  and  spits,  and  at  ephemeral 
river  mouths  and  coastal  lagoons  and  ponds  either  from  the  seaward  side  (Kraus 
et  al.  2002)  or  from  the  bay  side  (Kraus  et  al.  2008)  due  to  sustained  high  water 
level  and  wave  overtopping.  The  main  causes  of  breaching  are  overtopping  flow 
and  waves,  as  well  as  piping  and  seepage  flows.  Breaches  can  cause  loss  of 
human  life,  property,  infrastructure,  and  transportation  corridors;  endanger 
navigation  and  stability  of  adjacent  inlets;  and  degrade  the  environment  by  the 
exposure  of  bay  perimeter  to  sea  waves  and  by  the  change  in  bay  salinity. 

Coastal  barrier  and  dike  breach  modeling  has  attracted  more  and  more  attention 
in  recent  years.  Visser  (1998),  Kraus  and  Hayashi  (2005)  and  D’Eliso  (2007) 
developed  simplified  physically-based  coastal  barrier  or  dike  breaching  models, 
in  which  the  breach  cross-section  is  approximated  as  a  trapezoid  or  rectangle, 
the  overtopping  flow  is  estimated  using  the  broad-crested  weir  equation,  the 
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Keulegan  equation  or  the  Bernoulli  equation,  and  the  breach  morphological 
change  is  determined  using  simplified  sediment  transport  models.  Parkinson  and 
Stretch  (2007)  investigated  the  temporal  evolution  of  sand  barrier  breaching  in 
estuaries  and  coastal  lagoons  based  on  laboratory  experiments,  and  suggested  a 
scaling  method  of  breach  formation  time  and  peak  outflow,  considering  the 
effects  of  outflow  volume,  hydraulic  head,  and  barrier  breadth.  Roelvink  et  al. 
(2009)  developed  a  depth-averaged  2-D  model  of  coastal  barrier  breaching  in 
which  the  shallow  water  equations  with  wave  action  and  the  non-equilibrium 
sediment  transport  equations  are  solved  in  a  decoupled  way  using  a  finite 
difference  method. 

Wu  et  al.  (2009)  developed  a  depth-averaged  2-D  numerical  model  for 
simulating  the  morphodynamic  processes  of  earthen  embankment  breaching 
caused  by  overtopping  flow.  It  uses  a  finite-volume  shock-capture  scheme  based 
on  the  HLL  approximate  Riemann  solver  on  rectangular  mesh  to  handle  the 
breach  flow  that  is  usually  in  mixed  flow  regimes  and  with  discontinuities.  It 
adopts  a  coupled  solution  algorithm,  particularly  considering  the  effects  of 
strong  sediment  transport  and  rapid  bed  changes  on  the  flow  during  the 
breaching  process.  This  paper  extends  the  model  of  Wu  et  al.  to  simulate  the 
breaching  processes  of  coastal  barriers.  The  governing  equations,  numerical 
methods  and  validation  of  the  developed  model  are  presented  in  the  next 
sections. 

Governing  Equations 

In  comparison  with  dam  and  levee  breaching  in  inland  waters,  coastal  barrier 
and  levee  breaching  processes  are  much  more  complex.  The  breaching  can  occur 
from  either  side  of  the  barrier,  and  it  is  controlled  by  tidal  flow,  surge,  wind  and 
waves.  A  coastal  barrier  breach  from  the  lagoon  side  usually  opens  rapidly  and 
prevents  wave  penetration  because  of  the  ebb  flow.  An  ephemeral  ebb  delta  can 
also  form.  After  water  level  in  the  bay  drops  and  the  exiting  flow  is  reduced, 
waves  can  penetrate  the  breach  at  higher  tide,  causing  a  second  rapid  expansion 
of  its  width  (Bokuniewicz  et  al.  201 1).  A  breach  from  the  sea  side  is  oftencaused 
by  storm  surge  and  waves.  The  setup  of  water  level  by  strong  onshore-directed 
wind  and  surge  during  a  storm  and  the  run-up  of  waves  contribute  to  the 
inundation,  and  the  presence  of  waves  in  the  incipient  breach  increases  sediment 
mobilization  and  transport.  A  barrier  breach  may  close  naturally  by  the  sediment 
transported  from  adjacent  beaches  and  shores,  or  it  can  increase  in  size  and 
become  a  new  inlet  or  estuary  (Kraus  and  Wamsley  2003). 

To  simplify  the  model  formulation,  this  study  focuses  on  the  initial  development 
stage  under  supercritical  flow,  which  is  the  most  critical  stage  of  a  coastal 
breach  during  overtopping  flow.  In  this  stage,  the  wave  radiation  stress. 
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turbulent  stress,  and  Coriolis  force  are  negligible  in  comparison  with  the  inertia 
and  gravitational  forces.  Because  the  breach  opens  rapidly,  the  interaction  of 
flow,  sediment  transport  and  channel  bed  is  significant;  thus,  the  following 
generalized  2-D  shallow  water  equations  are  herein  used  (Wu  2008): 


d(pli)  d(puh)  d(pvh)  dzu 

~ir+—+^r+p^=0 


^-(puh)  +  ^-(pu2h)  +  ^-(puvh)  +  pgh^  +  ^-gh2^-  +  pg 
at  ox  ay  ox  1  ox  n 


n  mbuU  _ 


(2) 


-(^h)  +  -(pu^)  +  -{^h)  +  pgh^  +  \gh^  +  pg  /;1/3 


n  mbvU  =  Q 


(3) 


where  t  is  the  time,  x  and  y  are  the  longitudinal  and  lateral  coordinates,  h  is  the 
flow  depth,  u  and  v  are  the  flow  velocities  in  x  and  y  directions,  U=(u2+v2)l~, 
dzb/dt  denotes  the  rate  of  change  in  bed  surface  elevation,  Zb  is  the  bed  surface 
elevation  above  a  reference  datum.  zs  is  the  water  level,  n  is  the  Manning 
roughness  coefficient,  g  is  the  gravitational  acceleration,  p  is  the  density  of  the 
water  and  sediment  mixture  in  the  water  column  determined  by  p=pw(  1  -C,)+psC, 
with  C,  being  the  volumetric  concentration  of  total-load  sediment,  pb  is  the 
density  of  the  water  and  sediment  mixture  in  the  bed  surface  layer  determined 
by  Pb=PwP'm+Ps(l-p'm  )  with  p'm  being  the  porosity  of  the  surface-layer  bed 
material,  pw  and  ps  are  the  water  and  sediment  densities,  and 
mh=[  I  +(dziJox)2+(oziJoy)2] 1/2  describes  the  effect  of  bed  slope. 


Simulation  of  the  morphodynamic  processes  under  barrier  breach  flow  can  be 
achieved  by  computing  bed  load  and  suspended  load  separately  or  the  total  load 
directly  in  the  case  of  non-cohesive  sediments.  Here,  the  latter  approach  is 
taken,  so  that  the  sediment  transport  and  bed  deformation  are  determined  with 
(Wu  2008) 


3(/iC, )  d(huCt)  d(hvC') 
dt  d x  d  y 


(4) 


i'-p'n) 


dt 


- ) 


(5) 


where  q ,»  is  the  total-load  sediment  transport  capacity,  and  L  is  the  non¬ 
equilibrium  adaptation  length  of  sediment  transport.  Note  that  the  coefficient  mb 
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in  Eqs.  (4)  and  (5)  considers  the  lateral  erosion  along  the  side  slope.  L  is  given 
by 


L  =  max 


(6) 


where  cos  is  the  settling  velocity  of  sediment  particles,  Lh  is  the  adaptation 
length  of  bed  load,  and  a  is  the  adaptation  coefficient  of  suspended  load.  Wu 
(2008)  summarizes  available  methods  to  evaluate  Lb  and  a  in  the  case  of 
sediment  transport  under  common  flows;  however,  whether  these  methods  are 
valid  in  the  case  of  barrier  breaching  needs  to  be  investigated. 

The  sediment  transport  capacity  q ,»  is  determined  using  the  formula  of  Wu  et  al. 
(2000),  with  the  modification  suggested  by  Wu  (2008)  to  consider  the  effect  of 
gravity  on  sediment  transport  over  steep  slope.  In  addition,  to  consider  the  effect 
of  sediment  concentration,  the  settling  velocity  cos  used  in  Wu  et  al.  (2000) 
formula  is  determined  by  Richardson  and  Zaki’s  (1954)  method:  cos=cos<)(  I  -Ct)m, 
with  m~4.0  and  cos 0  being  the  settling  velocity  of  single  sediment  particles  in 
quiescent,  distilled  water. 

Numerical  Methods 

The  governing  equations  are  solved  with  an  explicit  finite  volume  method  on  a 
non-staggered  (collocated)  rectangular  grid  system.  The  HLL  approximate 
Riemann  solver  (Harten  et  al.  1983)  is  applied  to  determine  the  fluxes  at  cell 
faces.  The  water  surface  gradient  is  discretized  by  forward  differencing  if  the 
velocity  is  positive  and  by  backward  differencing  if  the  velocity  is  negative.  The 
numerical  method  has  first-order  accuracy  in  both  time  and  space.  The  sliding  of 
loose  materials  along  the  two  side  walls  of  the  breach  is  handled  by  the  non- 
cohesive  slope  avalanching  algorithm  presented  by  Wu  (2008).  The  flow  and 
sediment  transport  equations  are  solved  in  a  coupled  fashion  at  each  time  step. 
Details  of  the  numerical  algorithm  are  discussed  in  Wu  et  al.  (2009). 

Because  the  above  solution  procedure  is  explicit,  the  computational  time  step 
should  be  limited  by  numerical  stability  conditions,  such  as  the  Courant- 
Friedrichs-Lewy  (CFL)  condition  for  flow  computation  and  additional 
conditions  for  sediment  transport  and  bed  change  computations.  By  trial  and 
error,  it  was  found  that  to  obtain  a  stable  solution,  the  bed  change  at  each  time 
step  should  be  less  than  about  10%  of  the  local  flow  depth.  In  this  study,  the 
time  step  length  is  determined  by  setting  a  Courant  number  of  0.25  and  limiting 
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the  ratio  of  bed  change  to  local  flow  depth  to  be  less  than  0.1.  This  implies  that 
the  length  of  the  time  step  varies  during  the  simulation. 

Model  Test 

Case  1.  Dike-Break  Flow  on  Movable  Bed  in  a  Sudden-Expanded  Channel 

Experiments  of  dike  or  dam  break  flow  over  a  movable  bed  were  conducted  at 
the  Civil  Engineering  Laboratory  of  the  Universite  Catholique  de  Louvain 
(UCL),  Belgium.  The  case  selected  to  test  the  present  model  was  performed  in  a 
6  m  long  flume  with  a  non-symmetrical  sudden  expansion  from  0.25  m  to  0.5  m 
in  width,  located  4  m  downstream  of  the  flume  entrance  (Fig.  1).  The  breaking 
of  the  dike  was  simulated  by  the  rapid  downward  movement  of  a  thin  gate  in  the 
middle  of  the  flume.  The  initial  bed  of  the  entire  flume  was  covered  by  a  0.1  m 
thick  horizontal  layer  of  fully  saturated  sand,  which  was  uniform  coarse  sand 
with  a  median  diameter  of  1.82  mm  and  specific  gravity  of  2.68,  deposited  with 
a  porosity  of  0.47.  The  initial  water  depth  was  0.25  m  on  the  upstream  side  of 
the  gate  and  0  m  in  the  downstream  side  (Palumbo  et  al.  2008). 
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Fig.  1.  Sketch  of  the  UCL  dike-break  flow  experiment  setup  (dashed  lines  and  circles  denote 
locations  of  measurement  cross-sections  and  ultrasonic  gauges,  respectively). 


The  temporal  evolution  of  the  water  level  was  measured  by  eight  ultrasonic 
gauges  placed  downstream  of  the  gate,  as  sketched  in  Fig.  1.  The  gauges  Ul, 
U3,  U5,  and  U7  were  aligned  along  a  line  0.125  m  away  from  the  right  flume 
wall  (wall  at  bottom  of  Fig.  1)  and  located  at  0.75,  1.2,  1.45,  and  1.95  m 
downstream  of  the  gate,  respectively.  The  gauges  U2,  U4,  U6,  and  U8  were 
aligned  along  a  line  0.125  m  away  from  the  left  wall  (top  wall  in  Fig.  1)  and 
located  at  1.2,  1.45,  1.95,  and  2.45  m  downstream  of  the  gate,  respectively.  The 
morphological  changes  generated  by  the  dike -break  flow  were  measured  at  nine 
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cross-sections  (numbered  CS1  to  CS9)  located  every  5  cm  from  10  to  50  cm 
after  the  sudden  expansion,  at  the  end  of  the  experiment  by  a  laser  sheet  imaging 
technique. 

The  computational  domain  is  represented  by  a  mesh  consisting  of  560x100 
nodes  in  longitudinal  and  lateral  directions.  The  longitudinal  grid  spacing  is 
uniformly  set  as  0.005  m  near  the  gate,  but  varies  in  the  regions  away  from  the 
gate.  The  lateral  grid  spacing  is  uniform  and  has  a  value  of  0.005  m.  The 
sediment  adaptation  length  is  determined  by  Eq.  (6),  with  the  suspended  load 
adaptation  coefficient  a  equal  to  4.0  and  the  bed  load  Lb  set  as  0. 1  m.  The  time 
step  length  is  determined  by  setting  a  Courant  number  of  0.25  and  limiting  the 
ratio  of  bed  change  to  local  flow  depth  to  be  less  than  0.1.  The  Manning 
roughness  coefficient  is  given  as  0.025. 

Fig.  2  shows  the  calculated  total  bed  changes  near  the  sudden  expansion.  One 
can  see  that  erosion  occurs  near  the  corner  of  the  expansion,  and  the  eroded 
sediment  deposits  mostly  at  the  expansion  zone.  Figs.  3(a)-(e)  compare  the 
measured  and  calculated  final  bed  levels  at  five  cross-sections:  CS1,  CS3,  CS5, 
CS7  and  CS9.  At  cross-sections  CS1  and  CS3,  both  the  experiment  and 
numerical  model  present  erosion  at  the  center  of  the  flume  near  the  corner  of  the 
sudden  expansion,  but  the  model  underestimates  the  deposition  at  the  expansion 
zone  along  the  left  flume  wall  and  overestimates  the  erosion  along  the  right  wall. 
At  cross-sections  CS5,  CS7,  and  CS9,  the  model  reproduces  well  the  bed 
profiles  observed  in  the  experiment.  The  deposition  and  erosion  magnitudes  and 
distributions  are  predicted  with  good  accuracy. 

Figs.  4(a)-(e)  compare  the  measured  and  calculated  water  levels  at  gauges  Ul, 
U2,  U3,  U6  and  U7.  One  can  see  that  the  model  reproduces  well  the  general 
trend  of  the  dike-break  wave  propagation.  On  the  five  gauges,  the  root-mean- 
square  error  of  the  calculated  water  levels  is  0.0089  m,  and  the  correlation 
coefficient  (R2)  between  the  measured  and  calculated  water  levels  is  0.92. 
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Fig.  2.  Calculated  total  bed  changes  near  the  sudden  expansion. 
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Case  2.  Sea  Dike  Breaching  in  a  Tidal  Inlet 

The  developed  model  has  been  further  tested  by  comparison  with  the  Zwin’94 
field  experiment  measurements  of  sea  dike  breaching  conducted  at  the  Zwin 
Channel,  a  tidal  inlet  near  the  Dutch-Belgian  border.  The  channel  connects  the 
nature  reserve  “Het  Zwin”  with  the  North  Sea,  and  the  reserve  surface  area  is 
about  1.5  knr,  consisting  largely  of  marshes  and  gullies.  The  mean  tidal  prism  is 
about  350,000  m3.  A  dike  or  barrier  was  built  with  local  sand  across  the  Zwin 
Channel,  and  its  crest  was  3.3  m  above  the  mean  sea  level.  Inclination  of  the 
outer  slope  (seaward  side)  was  1:1.6  and  that  of  the  inner  slope  (bay  side)  1:3. 
The  width  of  the  crest  was  about  8.0  m.  The  bottom  elevation  of  the  Zwin 
Channel  near  its  axis  and  outside  the  sand  dike  varied  between  0.5  and  0.9  m, 
with  an  average  of  0.7  m  above  the  mean  sea  level.  The  dike  sediments 
consisted  of  ‘original’  Zwin  sand  with  a  median  diameter  of  0.185  mm  and 
additional  sand  with  a  median  diameter  of  0.315  mm  transported  alongshore 
from  the  Belgium  coast.  The  initial  breach  had  a  depth  of  0.8  m,  and  a  width  of 
about  1.0  m  near  its  bottom  or  about  3.6  m  near  the  crest  of  the  dike.  The  tidal 
level  outside  the  dike  was  measured  during  the  breaching,  and  the  highest  water 
level  was  0.27  m  above  the  initial  breach  bottom.  The  experiment  was  done  in 
quiet  autumn  weather  with  wind-speed  of  about  2  m/s  and  negligible  wave 
height  against  the  sand  dike  (Visser  1998). 

The  present  study  could  not  simulate  the  exact  breach  process  because  the 
topography  in  the  polder  area  is  not  clear.  Therefore,  the  simulation  considers  a 
hypothetical  case  with  a  sand  dike  with  the  same  geometry  and  material  as  used 
in  the  experiment.  The  channel  is  approximated  as  a  straight  channel  with  a 
width  of  60  m  (about  20  m  wider  than  the  maximum  breach  width).  The  sand 
dike  is  located  at  30  m  from  the  inlet.  The  channel  has  approximately  the  same 
storage  and  surface  area  of  the  polder  as  those  estimated  by  Smit  and  Snip 
(1995). 

A  mesh  consisting  of  500x200  nodes  in  longitudinal  and  lateral  directions  is 
used  in  the  simulation.  The  mesh  is  uniform  near  the  sand  dike,  but  nonuniform 
in  the  channel  away  from  the  dike.  The  grid  spacing  near  the  sand  dike  is  0.24  m 
in  the  longitudinal  direction  and  0.3  m  in  the  lateral  direction.  The  dike  material 
porosity  is  estimated  as  0.4.  The  sediment  adaptation  length  L  is  calibrated  as 
2.0  m.  A  representative  diameter  of  about  0.3  mm  is  specified  for  the  sand  in 
this  simulation.  The  bed  material  repose  angle  is  set  to  be  32°.  The  seaward  side 
of  the  channel  is  set  as  an  open  boundary,  with  the  water  levels  specified  as  the 
measured  ones  and  the  velocity  copied  from  the  adjacent  internal  nodes.  A  solid 
wall  boundary  condition  is  applied  at  the  other  end  and  along  two  banks  of  the 
channel. 
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Figs.  5(a)-(e)  show  the  simulated  flow  pattern  and  bed  topography  at  elapsed 
times  of  1.5,  4,  8.5,  16,  and  40  min  after  breaching,  respectively.  One  can  see 
the  evolution  of  the  breach  with  time.  As  the  flow  overtops  the  dike  crest  at  the 
initial  breach,  erosion  starts  at  the  downstream  (inner)  slope  and  develops 
upstream,  resulting  in  decrease  of  the  width  of  the  dike  crest.  As  the  bed  is 
eroded  vertically,  the  sediments  on  the  two  side  walls  of  the  breach  slide  at  the 
repose  angle  and  the  breach  is  widened  gradually.  A  large  scour  hole  appears  at 
the  center  of  the  channel  near  the  breach,  while  the  eroded  sediments  are 
transported  by  the  flow  and  immediately  deposited  along  the  sides  and 
downstream  edge  of  the  scour  hole.  The  numerical  model  gives  an  intuitively 
reasonable  general  flow  pattern  and  change  in  depth  contours  near  the  breach. 
The  model  produces  antidunes  or  standing  waves  on  the  bottom  at  the  early 
stage  of  breaching  with  supercritical  flow,  and  these  disappear  at  the  end  of 
breaching  after  the  flow  reduces  and  becomes  subcritical. 


Fig.  5.  Simulated  flow  pattern  and  breach  topography  (background  as  bed  elevation  contours  in 
meter;  x  and  y  coordinates  in  meter). 
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Fig.  6  compares  the  measured  and  simulated  breach  widths.  The  breach  width 
starts  to  increase  significantly  after  t=5.5  min  but  increases  slowly  after  about  35 
min.  The  breaching  process  finishes  at  about  53  min  when  the  water  level  in  the 
bay  rises  to  approximately  the  same  as  the  tidal  level  outside  the  sand  dike.  The 
measured  breach  width  at  60  min  is  41  m.  while  the  simulated  one  is  40.5  m. 
The  simulation  results  agree  generally  well  with  the  measurements 


Time  (min) 


Fig.  6.  Comparison  of  measured  and  simulated  breach  widths. 


Conclusions 

A  depth-averaged  2-D  finite-volume  model  was  established  to  simulate  the 
coastal  barrier  breaching  processes  accompanying  overtopping  flows.  The 
effects  of  sediment  concentration  and  bed  change  on  the  flow  are  represented  in 
the  continuity  and  momentum  equations.  The  sediment  transport  model 
simulates  the  non-equilibrium  transport  of  total  load.  An  explicit  algorithm 
based  on  the  finite  volume  method  is  adopted  to  solve  flow,  sediment  transport, 
and  bed  change  equations.  The  model  adopts  the  HLL  approximate  Riemann 
solver  to  describe  the  mixed-regime  flows  near  the  breach.  To  enhance  the 
performance  of  the  model,  a  varying  time  step  is  implemented  that  satisfies  the 
CFL  condition  and  ensures  the  bed  change  at  each  time  step  is  less  than  ten 
percent  of  the  local  flow  depth. 

The  developed  model  was  tested  using  a  laboratory  experiment  of  dike-break 
flow  on  a  movable  bed  carried  out  at  the  Universite  Catholique  de  Louvain,  and 
a  field  experiment  of  sea  dike  breaching  conducted  at  the  Zwin  Channel,  which 
is  a  tidal  inlet  near  the  Dutch-Belgian  border.  The  model  overall  predicts  well 
the  breach  width  and  the  flow  pattern  near  the  breach.  The  model  simulates  that 
antidunes  or  standing  waves  appear  in  the  early  stage  of  breaching  with 
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supercritical  flow,  and  they  disappear  at  the  end  of  breaching  when  the  flow 
becomes  subcritical.  Replication  of  these  processes  implies  that  the  model  is 
capable  of  simulating  the  initial  (catastropic)  developmental  stage  of  a  breach, 
when  its  depth  and  width  change  rapidly  under  supercritical  flow.  To  better 
represent  the  complex  breaching  of  barrier  and  dike  in  the  coastal  context,  the 
model  will  be  enhanced  in  the  near  future  to  take  into  account  the  effect  of 
waves. 
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